function [CRL,CRR,CRB,CRT,totLL,totLR,totLB,totLT] = crespL(phi,dx,dy,dc,mt,xfm,yfm,numg)
% current response generator

% number of fine meshes in x and y coordinates
N = sum(xfm); 
M = sum(yfm);
% number of coarse meshes in x and y coordinates
CN = length(xfm);
CM = length(yfm); poo=0.5;
LL = zeros(M+1,1); LR = LL; LB = zeros(N+1,1); LT = LB;
for g = 1:numg
    
% for L and R
for j = 1:M+1
    for ci = 1:CM
        if j-1 <= sum(yfm(1:ci)), mty = ci; break; end
    end
    for ci = 1:CM
        if j  <= sum(yfm(1:ci)), mtyu = ci; break; end
    end  
    % compute left response:  Jout = 0.25*phi + 0.5*D*(phi'-phi)/dx    
    i = 1;
    k = i+(j-1)*( N+1 ); ko = (i+1)+(j-1)*( N+1 );
    CRL( j, g) = 0.00*phi( k, g) + ...
                poo * ( dc(mt(1,mty),g) + dc(mt(1,mtyu),g) ) / dx(i) * ...
                ( phi( ko, g)-phi( k, g) );      
    LL(j) = LL(j) + 0.5*( dc(mt(1,mty),g) + dc(mt(1,mtyu),g) ) * ...
           (phi( ko, g)-phi( k, g))/dx(i);        
    % compute right response
    i = N+1;
    k = i+(j-1)*( N+1 ); ko = (i-1)+(j-1)*( N+1 );
    CRR( j, g) = 0.00*phi( k, g) + ...
                poo * ( dc(mt(end,mty),g) + dc(mt(end,mtyu),g) ) / dx(i-1) * ...
                ( phi( ko, g)-phi( k, g) );      
    LR(j) = LR(j) + 0.5*( dc(mt(end,mty),g) + dc(mt(end,mtyu),g) )* ...
           (phi( ko, g)-phi( k, g))/dx(i-1);          
end

% for B and T
for i = 1:N+1
    for ci = 1:CN
        if i-1 <= sum(xfm(1:ci)), mtx = ci; break; end
    end
    for ci = 1:CN
        if i  <= sum(xfm(1:ci)), mtxo = ci; break; end
    end
    % compute bottom response
    j = 1;
    k = i+(j-1)*( N+1 );  ko = (i-1)+(j)*( N+1 )+1;
    CRB( i, g) = 0.00*phi( k, g) + ...
                poo * ( dc(mt(mtx,1),g) + dc(mt(mtxo,1),g) ) / dy(j) * ...
                ( phi( ko, g)-phi( k, g) );      
    LB(i) = LB(i) + 0.5*( dc(mt(mtx,1),g) + dc(mt(mtxo,1),g) )* ...
            (phi( ko, g)-phi( k, g))/dy(j);          
    % compute top response        
    j = M+1;
    k = i+(j-1)*( N+1 ); ko = (i-1)+(j-2)*( N+1 )+1;
    CRT( i, g) = 0.00*phi( k, g) + ...
                poo * ( dc(mt(mtx,end),g) + dc(mt(mtxo,end),g) ) / dy(j-1) * ...
                ( phi( ko, g)-phi( k, g) );     
    LT(i) = LT(i) + 0.5*( dc(mt(mtx,end),g) + dc(mt(mtxo,end),g) )* ...
            (phi( ko, g)-phi( k, g))/dy(j-1);         
end

end
dxx = [0.5*dx(1) dx(1:end-1)' 0.5*dx(end)]';
dyy = [0.5*dy(1) dy(1:end-1)' 0.5*dy(end)]';
totLL = sum( dyy(2:end-1).*( LL(2:end-1) ) );
totLR = sum( dyy(2:end-1).*( LR(2:end-1) ) );
totLB = sum( dxx(2:end-1).*( LB(2:end-1) ) );
totLT = sum( dxx(2:end-1).*( LT(2:end-1) ) );

